function [u_x,u_y] = ex33_AnalyticalDisplacement(x,y,a,mu,k)
%function [u_x,u_y] = ex33_AnalyticalDisplacement(x,y,a,mu,k)
%
% x,y = nodal coordinates
% a   = radius of the hole
%

eps   = 1.0e-5;
a     = abs(a);
r     = sqrt(x*x+y*y);

% mu    = E/2/(1+nu);
% k     = 3-4*nu;        % plane strain
% %k     = (3-nu)/(1+nu); % plane stress

aux = a/8/mu;

if abs(r-a)<eps
    r = a;
end

if abs(y)<eps
    theta = 0;
elseif abs(x)<eps
    theta = 90*pi/180;
else    
    theta = atan(y/x);
end

u_x = aux*( r/a*(k+1)*cos(theta) + ...
            2*a/r*((1+k)*cos(theta) + cos(3*theta)) - ...
            2*a*a*a/r/r/r*cos(3*theta) );
          
u_y = aux*( r/a*(k-3)*sin(theta) + ...
            2*a/r*((1-k)*sin(theta) + sin(3*theta)) - ...
            2*a*a*a/r/r/r*sin(3*theta) );

